module testmoments
implicit none
contains
function ejm2(b,eta,si)
!E(jim*jih)
use normalmoms
use linear_operators
double precision eta,si
double precision ejm2,b(:)
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,p1,p2
integer ind

call paramstest(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
p1=h4-4*h3+6*h2-4*h1+1 !E((h-1)^4)
p2=h3-2*h2+h1 !E(h*(h-1)^2)
!!
ejm2=c**4.0*p1*dot_product(b,b)*dot_product(b,b) &
	&+c**2.0*p2*(dot_product(b,b)*quadnorm1(am .tx. am) +4*dot_product(b,matmul( (am .xt. am), b))&
	&+dot_product(b,b)*quadnorm1(am .tx. am))+h2*quadnorm2((am .tx. am), (am .tx. am))
end function ejm2
!!
function ejmem(b,eta,si)
!E(jih*em)
use normalmoms
use linear_operators
double precision eta,si
double precision b(:)
double precision ejmem(1:size(b))
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,p1,p2
integer ind

call paramstest(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
p1=h3-1-3*h2+3*h1 !E((h-1)^3)
p2=h2-h1 !E(h*(h-1))
!!
ejmem=c**3.0*p1*dot_product(b,b)*b&
	&+2.0*c*p2*matmul(am.xt.am,b)&
	&+c*p2*quadnorm1(am .tx. am)*b
end function ejmem
!!
function ejmememp(b,eta,si)
!E(jih*em*em')
use normalmoms
use linear_operators
use ghs
double precision eta,si
double precision b(:)
double precision ejmememp(1:size(b),1:size(b))
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,p1,p2
integer ind

call paramstest(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
p1=h4-4*h3+6*h2-4*h1+1 !E((h-1)^4)
p2=h3-2*h2+h1 !E(h*(h-1)^2)

ejmememp=c**4.0*p1*dot_product(b,b)*matmul(matvec(b),.t. matvec(b))&
	   &+c**2.0*p2*dot_product(b,b)*(am .xt. am)&
	   &+2.0*c**2.0*p2*((matmul(matvec(b),.t.matvec(b)).x.am) .xt. am)&
	   &+2.0*c**2.0*p2*((am.xt.am) .x. matmul(matvec(b),.t.matvec(b)))&
	   &+c**2.0*p2*matmul(matvec(b),.t.matvec(b))*quadnorm1(am .tx. am)&
	   &+h2*((am .x. quadvec1(am .tx. am)).xt. am)

end function ejmememp
!!
function ejm2ememp(b,eta,si)
!E(jih*em*em')
use normalmoms
use linear_operators
use ghs
double precision eta,si
double precision b(:)
double precision ejm2ememp(1:size(b),1:size(b))
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,h5,h6,p1,p2,p3
integer ind

call paramstest(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
h5=eh(nu,gam,5)
h6=eh(nu,gam,6)
p1=h6-6*h5+15*h4-20*h3+15*h2-6*h1+1 !E((h-1)^6)
p2=h5-4*h4+6*h3-4*h2+h1 !E(h(h-1)^4)
p3=h4-2*h3+h2 !E(h^2(h-1)^2)

ejm2ememp=c**6.0*p1*(dot_product(b,b)**2.0)*matmul(matvec(b),.t.matvec(b))&
		&+4.0*c**4.0*p2*matmul(matvec(b),.t.matvec(b))*dot_product(matmul(.t. am,b),matmul(.t. am,b))&
		&+c**2.0*p3*matmul(matvec(b),.t.matvec(b))*quadnorm2(am.tx.am,am.tx.am)&
		&+2.0*c**4.0*p2*dot_product(b,b)*matmul(matvec(b),.t.matvec(b))*quadnorm1(am .tx. am)&
		&+4.0*c**4.0*p2*dot_product(b,b)*((matvec(b).xt.matvec(b)).x. (am .xt. am))&
		&+4.0*c**2.0*p3*((((matvec(b).xt.matvec(b)).x. am).x. quadvec1(am .tx. am)).xt. am)&
		&+4.0*c**4.0*p2*dot_product(b,b)*((am .xt. am) .x. (matvec(b) .xt. matvec(b)))&
		&+4.0*c**2.0*p3*(((am .x. quadvec1(am .tx. am)) .xt. am) .x. (matvec(b) .xt. matvec(b)))&
		&+c**4.0*p2*(dot_product(b,b)**2.0)*(am .xt. am)&
		&+4.0*c**2.0*p3*((am .x. quadvec1((am .tx. matvec(b)).xt.(am .tx. matvec(b)))) .xt. am)&
		&+h3*((am .x. quadvec2(am .tx. am, am.tx.am)).xt.am)&
		&+2*c**2.0*p3*dot_product(b,b)*((am .x. quadvec1(am .tx. am)).xt. am)

end function ejm2ememp
!!
function ejm2em(b,eta,si)
!E(jih^2*em)
use normalmoms
use linear_operators
use ghs
double precision eta,si
double precision b(:)
double precision ejm2em(1:size(b))
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,h5,p1,p2,p3
integer ind

call paramstest(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
h5=eh(nu,gam,5)
p1=h5-5*h4+10*h3-10*h2+5*h1-1 !E((h-1)^5)
p2=h4-3*h3+3*h2-h1 !E(h*(h-1)^3)
p3=h3-h2 !E(h^2*(h-1))

ejm2em=c**5.0*p1*(dot_product(b,b)**2.0)*b&
    &+4.0*c**3.0*p2*dot_product(matmul(.t.am,b),matmul(.t.am,b))*b&
	&+c*p3*quadnorm2(am.tx.am,am.tx.am)*b&
	&+2.0*c**3.0*p2*dot_product(b,b)*quadnorm1(am .tx.am)*b&
	&+4.0*c**3.0*p2*dot_product(b,b)*matmul(am .xt. am,b)&
	&+4.0*c*p3*matmul((am .x. (quadvec1(am .tx. am) .xt. am )),b)

end function ejm2em
!!
function ejm3em(b,eta,si)
!E(jih^2*jim*em)
use normalmoms
use linear_operators
use ghs
double precision eta,si
double precision b(:)
double precision ejm3em(1:size(b))
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,h5,h6,h7,p1,p2,p3,p4
integer ind

call paramstest(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
h5=eh(nu,gam,5)
h6=eh(nu,gam,6)
h7=eh(nu,gam,7)

p1=h7-7*h6+21*h5-35*h4+35*h3-21*h2+7*h1-1 !E((h-1)^7)
p2=h6-5*h5+10*h4-10*h3+5*h2-h1 !E(h(h-1)^5)
p3=h5-3*h4+3*h3-h2 !E(h^2(h-1)^3)
p4=h4-h3 !E(h^3(h-1))

ejm3em=c**7.0*p1*dot_product(b,b)**2.0*dot_product(b,b)*b&
       &+c**5.0*p2*dot_product(b,b)*dot_product(b,b)*quadnorm1(am.tx.am)*b&
	   &+2.0*c**5.0*p2*dot_product(b,b)*dot_product(b,b)*matmul(am.xt.am,b)&
	   &+4.0*c**5.0*p2*dot_product(b,b)*dot_product(matmul(.t.am,b),matmul(.t.am,b))*b&
	   &+2.0*c**5.0*p2*dot_product(b,b)**2.0*matmul(am.xt.am,b)&
	   &+2.0*c**3.0*p3*dot_product(b,b)*matmul((am.x.quadvec1(am.tx.am)).xt.am,b)& 
	   &+c**5.0*p2*dot_product(b,b)**2.0*quadnorm1(am.tx.am)*b&
	   &+c**3.0*p3*dot_product(b,b)*quadnorm2(am.tx.am,am.tx.am)*b& 
	   &+2.0*c**3.0*p3*dot_product(b,b)*matmul((am.x.quadvec1(am.tx.am)).xt.am,b)&
	   &+4.0*c**5.0*p2*dot_product(b,b)*dot_product(matmul(.t.am,b),matmul(.t.am,b))*b&
	   &+2.0*c**5.0*p2*dot_product(b,b)*dot_product(b,b)*matmul(am.xt.am,b)&
	   &+2.0*c**3.0*p3*dot_product(b,b)*matmul((am.x.quadvec1(am.tx.am)).xt.am,b)&
	   &+4.0*c**5.0*p2*dot_product(b,b)*dot_product(matmul(.t.am,b),matmul(.t.am,b))*b&
	   &+4.0*c**3.0*p3*quadnorm2(am.tx.am,(am.tx.(matvec(b).xt.matvec(b))).x.am)*b&          
	   &+8.0*c**3.0*p3*matmul((am.x.quadvec1((am.tx.(matvec(b).xt.matvec(b))).x.am)).xt.am,b)&
	   &+4.0*c**3.0*p3*quadnorm2(am.tx.am,(am.tx.(matvec(b).xt.matvec(b))).x.am)*b&
	   &+2.0*c**3.0*p3*dot_product(b,b)*matmul((am.x.quadvec1(am.tx.am)).xt.am,b)& !ok
	   &+2.0*c*p4*matmul((am.x.quadvec2(am.tx.am,am.tx.am)).xt.am,b)&
	   &+c**5.0*p2*dot_product(b,b)*dot_product(b,b)*quadnorm1(am.tx.am)*b&
	   &+c**3.0*p3*dot_product(b,b)*quadnorm2(am.tx.am,am.tx.am)*b&
	   &+2.0*c**3.0*p3*dot_product(b,b)*matmul((am.x.quadvec1(am.tx.am)).xt.am,b)&
	   &+4.0*c**3.0*p3*quadnorm2(am.tx.am,(am.tx.(matvec(b).xt.matvec(b))).x.am)*b&
	   &+2.0*c**3.0*p3*dot_product(b,b)*matmul((am.x.quadvec1(am.tx.am)).xt.am,b)&
	   &+2.0*c*p4*matmul((am.x.quadvec2(am.tx.am,am.tx.am)).xt.am,b)&
	   &+c**3.0*p3*dot_product(b,b)*quadnorm2(am.tx.am,am.tx.am)*b&
	   &+c*p4*quadnorm3(am.tx.am,am.tx.am,am.tx.am)*b&
	   &+2.0*c*p4*matmul((am.x.quadvec2(am.tx.am,am.tx.am)).xt.am,b)
end function ejm3em


function ejm3(b,eta,si)
!E(jih*jim^2)
use normalmoms
use linear_operators
use ghs
double precision eta,si
double precision b(:)
double precision ejm3
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,h5,h6,p1,p2,p3
integer ind

call paramstest(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
h5=eh(nu,gam,5)
h6=eh(nu,gam,6)
p1=h6-6*h5+15*h4-20*h3+15*h2-6*h1+1 !E((h-1)^6)
p2=h5-4*h4+6*h3-4*h2+h1 !E(h(h-1)^4)
p3=h4-2*h3+h2 !E(h^2(h-1)^2)

ejm3=c**6.0*p1*dot_product(b,b)**2.0*dot_product(b,b)&
     &+c**4.0*p2*dot_product(b,b)**2.0*quadnorm1(am.tx.am)&
	 &+8.0*c**4.0*p2*dot_product(b,b)*quadnorm1((am.tx.(matvec(b).xt.matvec(b))).x.am)&
	 &+2.0*c**4.0*p2*dot_product(b,b)*dot_product(b,b)*quadnorm1(am.tx.am)&
	 &+2.0*c**2.0*p3*dot_product(b,b)*quadnorm2(am.tx.am,am.tx.am)&
	 &+4.0*c**4.0*p2*dot_product(b,b)*quadnorm1((am.tx.(matvec(b).xt.matvec(b))).x.am)&
	 &+4.0*c**2.0*p3*quadnorm2((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am)&
	 &+8.0*c**2.0*p3*quadnorm2((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am)&
	 &+c**2.0*p3*dot_product(b,b)*quadnorm2(am.tx.am,am.tx.am)&
	 &+h3*quadnorm3(am.tx.am,am.tx.am,am.tx.am)
end function ejm3

function ejm4(b,eta,si)
!E(jih^2*jim^2)
use normalmoms
use linear_operators
use ghs
double precision eta,si
double precision b(:)
double precision ejm4
double precision nu,gam,c,am(1:size(b),1:size(b))
double precision h1,h2,h3,h4,h5,h6,h7,h8,p1,p2,p3,p4
integer ind

call paramstest(b,am,c,nu,gam,eta,si)
!!
h1=1.0d0
h2=eh(nu,gam,2)
h3=eh(nu,gam,3)
h4=eh(nu,gam,4)
h5=eh(nu,gam,5)
h6=eh(nu,gam,6)
h7=eh(nu,gam,7)
h8=eh(nu,gam,8)

p1=h7-6*h6+15*h5-20*h4+15*h3-6*h2+h1 !E(h(h-1)^6)
p2=h6-4*h5+6*h4-4*h3+h2  !E(h^2(h-1)^4)
p3=h8-8*h7+28*h6-56*h5+70*h4-56*h3+28*h2-8*h1+1   !E((h-1)^8)
p4=h5-2*h4+h3  !E(h^3(h-1)^2)

ejm4=16.0*c**6.0*p1*dot_product(b,b)*dot_product(b,b)*quadnorm1((am.tx.(matvec(b).xt.matvec(b))).x.am)&
      &+c**4.0*p2*dot_product(b,b)**2.0*quadnorm2(am.tx.am,am.tx.am)&
	  &+4.0*c**6.0*p1*dot_product(b,b)**2.0*quadnorm1((am.tx.(matvec(b).xt.matvec(b))).x.am)&
	  &+h4*quadnorm4(am.tx.am,am.tx.am,am.tx.am,am.tx.am)&
	  &+2.0*c**6.0*p1*dot_product(b,b)**2.0*dot_product(b,b)*quadnorm1(am.tx.am)&
	  &+2.0*c**6.0*p1*dot_product(b,b)*dot_product(b,b)**2.0*quadnorm1(am.tx.am)&
	  &+4.0*c**6.0*p1*dot_product(b,b)**2.0*quadnorm1((am.tx.(matvec(b).xt.matvec(b))).x.am)&
	  &+c**8.0*p3*dot_product(b,b)**2.0*dot_product(b,b)**2.0&
	  &+16.0*c**4.0*p2*dot_product(b,b)*quadnorm2((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am)&
	  &+8.0*c**4.0*p2*dot_product(b,b)*quadnorm2((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am)&
	  &+16.0*c**4.0*p2*quadnorm2((am.tx.(matvec(b).xt.matvec(b))).x.am,(am.tx.(matvec(b).xt.matvec(b))).x.am)&
	  &+2.0*c**2.0*p4*dot_product(b,b)*quadnorm3(am.tx.am,am.tx.am,am.tx.am)&
	  &+4.0*c**4.0*p2*dot_product(b,b)*dot_product(b,b)*quadnorm2(am.tx.am,am.tx.am)&
	  &+8.0*c**4.0*p2*dot_product(b,b)*quadnorm2((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am)&
	  &+2.0*c**2.0*p4*dot_product(b,b)*quadnorm3(am.tx.am,am.tx.am,am.tx.am)&
	  &+4.0*c**2.0*p4*quadnorm3((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am,am.tx.am)&
	  &+c**4.0*p2*dot_product(b,b)**2.0*quadnorm2(am.tx.am,am.tx.am)&
	  &+16.0*c**4.0*p2*dot_product(b,b)*quadnorm2((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am)&
	  &+4.0*c**2.0*p4*quadnorm3((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am,am.tx.am)&
	  &+16.0*c**2.0*p4*quadnorm3((am.tx.(matvec(b).xt.matvec(b))).x.am,am.tx.am,am.tx.am)
end function ejm4
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine paramstest(b,am,c,nu,gam,eta,si)
use linear_operators
use bessels
use ghs
double precision, intent(in):: b(:),eta,si
double precision, intent(out):: nu,gam,c,am(1:size(b),1:size(b))
double precision a,normb
integer k,n
!ind=1 normal
!ind=2 inverted
!ind=3 all b
!ind=4 all b

nu=-1.0d0/(2.0d0*eta)
gam=(1-si)/si
c=ghcfun0(b,eta,si)
if (dot_product(b,b).gt.1.0d-3) then
	am=eye(size(b))+((c-1)/dot_product(b,b))*(matvec(b) .xt.matvec(b))
	am=.t.chol(am)
else 
	a=besseld(nu+1.0d0,gam)-1.0d0
	normb=dot_product(b,b)
	am=eye(size(b))+(-a+2.0d0*(a*normb)**2.0d0-5.0d0*a**3.0d0*normb**2.0d0)*(matvec(b) .xt. matvec(b))
	am=.t.chol(am)
end if

end subroutine paramstest

function eh(nu,gam,k)
!E(ht**k)
use bessels
double precision nu,gam,eh
integer k,i
eh=1
do i=1,k
	eh=eh*besselr(nu+i-1,gam)/besselr(nu,gam)
end do

end function eh
end module testmoments